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We explore the Renyi entanglement entropies of a one-dimensional (line) subsystem of length L 
embedded in two-dimensional L x L square lattice for quantum spin models whose ground-state 
breaks a continuous symmetry in the thermodynamic limit. Using quantum Monte Carlo simula¬ 
tions, we first study the Ji — J 2 Heisenberg model with antiferromagnetic nearest-neighbor Ji > 0 
and ferromagnetic second-neighbor couplings J 2 < 0. The signature of SU(2) symmetry breaking 
on finite size systems, ranging from L = 4 up to L = 40 clearly appears as a universal additive 
logarithmic correction to the Renyi entanglement entropies: Z^lnL with l q ~ 1, independent of the 
Renyi index and values of J 2 . We confirm this result using a high precision spin-wave analysis (with 
restored spin rotational symmetry) on finite lattices up to 10 s x 10 5 sites, allowing to explore further 
non-universal finite size corrections and study in addition the case of U(l) symmetry breaking. Our 
results fully agree with the prediction l q = ug/2 where ug is the number of Goldstone modes, by 
Metlitski and Grover [arXiv: 1112.5166]. 


I. INTRODUCTION 

Entanglement entropy (EE) is now well recognized as a 
very powerful tool to diagnose various quantum states of 
matter 1,2 . For interacting quantum systems in dimension 
D > 2, the ground-state EE of a given spatial partition 
A embedded in a larger system scales with the perimeter 
lA of A, following the so-called area-law 3,4 for any Renyi 
index q > 0 

S q = Y^^hi^TrfpA] 9 ) = a q i A A - (1-1) 

where pa is the reduced density matrix of the subsystem 
A. While the leading part a q (-A is not expected to reflect 
the universality of the phase, sub-leading terms (the ellip¬ 
sis in Eq. (1.1) above) may encode it, as first discovered 
for topological order 5,6 . For systems which exhibit a true 
long-range order in the ground-state with a continuous 
symmetry breaking in the thermodynamic limit, the EE 
of a subsystem has been predicted' to exhibit a univer¬ 
sal additive logarithmic correction to the area law term, 
with a prefactor in dimension D = 2 controlled by the 
number of Goldstone modes hq associated to the broken 
symmetry: 


Sq — dq^A H—bit/4 + • • ■ (1-2) 



Figure 1. Schematic picture for the Ji — J 2 square lattice. 
The line shaped subsystem A of length £a is shown in green. 


not agree with the prediction. Also, the expected fac¬ 
tor of two between the logarithmic terms for SU(2) and 
U(l) was not clearly observed 11 . The possible reasons 
for such discrepancies are temperature and/or statisti¬ 
cal effects, as well as the importance of further finite-size 
corrections (beyond the log term) which might be hard 
to capture with QMC simulations on finite-size systems. 
Also one has to carefully subtract contributions to loga¬ 
rithmic corrections to S q from corners if present in the 
subsytem geometry 10,12 : these contributions are usually 
numerically small and their estimates from QMC simula¬ 
tions are suffering from the above mentioned difficulties. 


In such finite systems (with N = L x L sites), there are 
two types of excitations: the Anderson tower of states 
(TOS) 8 with an energy scaling as 1/L 2 and no Goldstone 
modes (SW for quantum magnets) with a linear disper¬ 
sion ~ 1/L, both being key contributions for the expected 
logarithmic corrections in Eq. (1.2)'. As first detected 
using a modified spin-wave (SW) approach for the SU(2) 
symmetric Heisenberg antiferromagnet on a square lat¬ 
tice 9 , subsequent quantum Monte Carlo (QMC) calcula¬ 
tions 10-12 have also been able to capture additive loga¬ 
rithmic corrections, while estimates of the prefactor did 


Very recently, Kulchytskyy et al. used an improved 
estimator for S 2 (for a half-torus subsystem) with QMC 
simulations of the spin-); XY model on the square lat¬ 
tice 13 which allowed them to get a quite precise estimate 
for the prefactor of the log correction ~ 0.5, fully con¬ 
sistent with no = 1 Goldstone boson associated to the 
breaking of U(l) symmetry. Nevertheless, to the best of 
our knowledge there is no numerical study demonstrating 
the universality of Eq. (1.2), such as its independence on 
the Renyi index q , details of microscopic Hamiltonian or 
type of continuous symmetry breaking. 
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Figure 2. (Color online) QMC results for the entanglement Renyi entropies of the Ji — J 2 Heisenberg model for J 2 = —1 (left 
two panels) and J 2 = —3 (right two panels). We show the prefactor of the logarithmic scaling term obtained by fits to the form 
S q = a q L + l q lnL + b q + c q /L over fit ranges [L m m, imax] as a function of L m i n , with L m ax = 40 for q = 2 and L max = 36 for 
q = 3,4. Our results are consistent with l q = 1 independent of J 2 and q. For J 2 = —3, we also show the EE S 2 lrect obtained 
by a direct mixed ensemble calculation using the method of Ref. 11. 


In this paper, we aim at going further to test the pre¬ 
diction Eq. (1.2) for various values of q and for two quan¬ 
tum spin models having different symmetries, using a 
one-dimensional ring of length £ A = L as subsystem A 
(see Fig. 1) embedded in a L x L torus. This is the sim¬ 
plest possible corner-free bipartition scaling with L where 
the universal logarithmic correction proportional to the 
number of Goldstone modes should be present. 

We explore the Renyi EEs S q for such a subsystem 
using two techniques: exact QMC simulations of S q with 
q = 2, 3,4 (Sec. II), and a semi-classical SW theory for 
finite size systems 14 ’ 15 where spin rotational invariance 
is restored such that both TOS and Goldstone modes 
are included (Sec. III). The choice of a line subsystem is 
advantageous for these two techniques: in QMC, we use 
the improved estimator introduced in Ref. 16, which is 
particularly efficient when subsystem volume is as small 
as possible (it is in fact minimal for the line subsystem), 
while the SW calculations are particularly simplified by 
full translation symmetry of the line subsystem, allowing 
for an analytical understanding of the In L term of 
Eq. (1.2). 

II. QUANTUM MONTE CARLO RESULTS 

For our quantum Monte Carlo calculations, we con¬ 
sider the spin-1/2 J\ — J 2 antiferromagnet defined on a 
bipartite L x L square lattice by the following Hamilto¬ 
nian 

nj.-j, = Sj + X! & • (2- 1 ) 

(ij) «b» 

where S are spin-| operators, interactions act between 
nearest neighbors (ij) and second nearest neighbours 


((ij)) along the diagonals of the square lattice (see 
Fig. 1). For this work, we consider antiferromagnetic 
nearest neighbors interactions J\ > 0 and ferromagnetic 
second neighbors interactions J 2 < 0, for which it is 
known that the ground-state exhibits antiferromagnetic 
long-range order, thus breaking SU(2) symmetry associ¬ 
ated with two Goldstone modes (independent of J 2 < 0). 
The motivation for adding the second neighbors interac¬ 
tion J 2 is to check the universality of the results with 
respect to microscopical variations of the Hamiltonian 
(different values of J 2 ) without changing the nature of 
the ground-state and of the low-lying excitations. Addi¬ 
tionally, as | J 2 | is increased, the antiferromagnetic long- 
range order is enhanced (i.e. larger values of the order 
parameter), and we therefore expect lower EEs as we get 
closer to a classical Heisenberg antiferromagnet. 

We perform extensive QMC simulations of this model 
for two different values of J 2 using the stochastic se¬ 
ries expansion (SSE) algorithm 1 '’ 18 . We compute the 
Renyi EE S q for q = 2,3,4 using a recently introduced 
decomposition 16 that benefits from subsystem symme¬ 
tries. This method is particularly useful when the surface 
of the subsystem scales as its volume, i.e. if the subsys¬ 
tem volume is minimal without introducing geometrical 
effects, such as corners. In this sense, this method is 
optimal for the line shaped subsystem in Fig. 1. 

Our simulations are performed in the finite tempera¬ 
ture formulation of the SSE at low enough temperatures 
in order to capture only ground-state physics. While the 
finite size gap of the tower of states scales as the inverse 
of the total number of spins N = L x L in the system, 
one would expect that it is necessary to scale the inverse 
temperature f3 linearly in N. For the system sizes we 
studied (up to L = 40), we find however that the re¬ 
sults for the EE of simulations at inverse temperatures 
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/3 = 8L and /3 = 4L agree within errorbars: we there¬ 
fore performed all calculations at inverse temperatures 
higher than /3 = 4 L. Fig. 2 shows the QMC result of the 
line EEs as a function of system size for different Renyi 
indices and two values of J 2 ■ For J 2 = — 3 and q = 2, 
we also perform an independent set of QMC simulations 
in an extended ensemble where EE is directly computed 
from the ratio of partition functions 11 . We obtain a per¬ 
fect agreement between the two methods. Note, that the 
decomposition of the EE as described in Ref. 16 allows 
us to access larger system sizes (in particular for higher 
Renyi indices) with very high precision. 

We fit our results for the line EEs to the scaling ansatz 

Sq = dqt lq\ll£ ~\~ bq Cq /E (2-2) 

to infer if the prefactor of the logarithmic term is in¬ 
deed l q = net 2 (= 1 for the ground-state of the model 
Eq. 2.1). We systematically reduce the fitting range 
[E m in, E max ] of included system sizes, always including 
the largest systems with E max ~ 40 and studying the 
best fit value of l q as a function of E m m as shown in the 
right panels of Fig. 2. Note that the errorbars stem from 
a careful bootstrap study of the stability of the fit intro¬ 
ducing gaussian resampling of the data and perturbations 
of the initial parameters. We have studied the statisti¬ 
cal behavior of the fit-data distance quantified by y 2 and 
find that the qualities 19 Q of the best fit are already very 
good (around 0.7... 0.9) using the scaling ansatz from 
equation (2.2). Hence, the quality of our data does not 
allow for an inclusion of higher order terms (which could 
result in overfitting statistical noise). 

While the log term is clearly present as nicely visible 
in the concavity of S q (L), we found it difficult to get a 
very precise estimate for the prefactor l q in the thermo¬ 
dynamic limit. What is clear however is that whereas 
the area law term does depend on the Renyi index q 
and J 2 , there is apparently no (^-dependence for l q . Tak¬ 
ing into account the largest E m i n , we can estimate that 
lq = 1.0(3), fully compatible with the prediction l q = 1, 
albeit with admittedly large error bars. In order to reach 
much larger systems (which would be helpful in studying 
the convergence of l q with E m in), we now consider a SW 
calculation of EE for the same setup of a line subsystem. 


III. SPIN-WAVE THEORY 

Modified SW theory for finite size systems 14,15 has 
been shown to be very useful for computing EEs of the 
square lattice Heisenberg antiferromagnet in Ref. 9. The 
crucial point is to artificially restore the spin rotational 
invariance in order to mimic the symmetric ground-state 
of a finite size system. For this, a size-dependent regular¬ 
izing external held h* is imposed to the system such that 
the SW-corrected order parameter is identically zero. 

We study the J\ — J 2 Heisenberg antiferromagnet 
Eq. (2.1) where SU(2) symmetry is restored by adding 


a small staggered magnetic held 

n . h - j2 = Sj + J2^2Si- + h* ]T(-i ys*, 

(ij) <(y» » 

(3.1) 

such that (Sf) = 0, as well as the ferromagnetic XY 
model 

hxy = -jJ 2 ( s i s J + s i s i)+ h *J2 s ?> ( 3 - 2 ) 

(ij) i 


where the transverse held h* is chosen such that (Sf) = 0 
and (Sf) = 0 in order to artificially restore the U(l) sym¬ 
metry. In both cases, the external held h* oc L~ 4 leads 
to a finite size gap A ~ y/h* ~ 1/E 2 , below the no lin¬ 
early dispersing Goldstone modes (no = 2 for the SU(2) 
Jl — J 2 model, and no = 1 for the U(l) XY model). This 
additional energy scale reproduces the TOS structure on 
hnite systems. For the analytical expressions below, we 
do not specify the value of the spin S, while for the nu¬ 
merical computations we explicitly consider S = 1/2. 

The calculation of the EEs in the SW approximation 
is eased by the quadratic nature of the SW Hamiltonian 
(at the linear harmonic level), as the reduced density ma¬ 
trix can be expressed as an exponential of a correlation 
matrix C involving only expectation values of two-point 
correlation functions 20,21 . The EEs of a subsystem com¬ 
posed of Na sites are obtained as a function of the Na 
eigenvalues i' 2 of this correlation matrix: 



The eigenvalues z/ 2 are typically obtained from numer¬ 
ical diagonalization of the matrix C, once its elements 
have been evaluated for the value h* (which depends on 
system size and parameters such as J 2 or the spin value 
S). In the case of a line subsystem, the numerical di¬ 
agonalization step can be circumvented by noticing that 
the translation invariance along the line implies that the 
matrix C is circulant 22 : its eigenvalues are given by the 
Fourier transform of its first line. Plugging in all the ex¬ 
pectation values and exploiting the convolution theorem, 
we obtain 


= J_ [ y- 2l(p,fc„) \ 

2 e ^ yjrf n(p, k y ) J 



B{p, k v ) 

Q(p, ky) 


2 

, (3.5) 


where p = —ir + ^-(j — 1) with j £ [1,E], and the SW 
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excitation spectrum f2(k) = yjA( k) 2 — B( k) 2 is given by 

h* 

A(k x , k y ) = 25J 2 cos k x cos + 2S(Ji — J 2 ) + — 
B(k X: k y ) = -SJi [cos + cos fc y ] (3.6) 

for the Ji — J 2 model, while for the XY case 

j h* 

A(k x , k y ) = —S'— [cosfca, + cos/c y ] + SJ + — 

B(k x ,k y ) = S— [cos/c^ +cosfc y ]. (3.7) 

The symmetry of the line subsystem allows to di¬ 
rectly compute EEs within the SW approximation us¬ 
ing Eq. (3.5) for systems of very large linear size (up 
to L ~ 10 5 ), which can hardly be reached if a nu¬ 
merical diagonalization of C is involved. Note that for 
such large systems, the regularizing field h* becomes ex¬ 
tremely small, and we resorted to arbitrary-precision nu¬ 
merics to ensure convergence of h* and corresponding v p . 

Considering now the case S = 1/2, our numerical re¬ 
sults for S q for different Renyi indices q are displayed in 
Fig. 3 (left) for the J\ — J 2 Heisenberg model (for different 
values of the diagonal coupling J 2 ) and in Fig. 4 (left) for 
the XY model. The precise value of S q being dominated 
by a non-universal area law term, one cannot directly 
compare the actual estimates of S q obtained within SW 
to our exact QMC results at small sizes due to the ap¬ 
proximations inherent to the SW approach. However, we 
expect the universal subleading logarithmic scaling term 
to be well captured by the modified SW theory. Indeed, 
fitting our SW data to the previous form Eq. (2.2) clearly 
yields an additive logarithmic term, as shown in Fig. 5 for 
the Ji — J 2 antiferromagnet with J 2 = — 1 and in Fig. 6 




Figure 3. (Color online) Left panel: EE S q of the line shaped 
subsystem in the Ji — J 2 model for different values of J 2 and 
Renyi indices q as obtained from the modified spinwave anal¬ 
ysis. Right panel: Prefactors l q of the logarithmic corrections 
obtained from fits of the form albl'l"d (c/. Eq. (3.8) for def¬ 
initions of the terms) as a function of the minimal size L m i n 
included in the fit. 




Figure 4. (Color online) Left panel: EE S q of the line shaped 
subsystem in the XY model for different Renyi indices q. 
Right panel: Logarithmic term in the scaling of the EE of 
a line in the XY model for different fit ranges [L m i n ,7 • 10 4 ] 
and different Renyi indices q, as obtained from a fit of the 
form albl'l"c!c (see Eq. (3.8)). 

for the XY model. The slow convergence of the coef¬ 
ficient of the logarithmic term suggests that subleading 
corrections beyond the log term in Eq. (2.2) have to be 
included. As we are not aware of any prediction for such 
subleading corrections, we perform fits using the general 
ansatz 

Sq = CLqL + Iqln L + b q 

+ l 2 In In L + l 3 In In In L + + c 1 , (3.8) 

Lj 1j 

leaving out systematically various terms. We use the 
shorthand notation a/M 2 / 3 c 1 c to label the various fit func¬ 
tions in the following figures (terms whose parameters do 
not appear in this string are not included in the fits). We 
find nonvanishing contributions for all terms and compar¬ 
ing carefully the distance of the fit to the data quantified 
by y 2 j it seems that the inclusion of all these terms yields 
the best fits. We show a representative analysis of dif¬ 
ferent fit functions in Fig. 5 for the Ji — J 2 model and 
in Fig. 6 for the XY model. The comparison of the dis¬ 
tances of the studied fit functions to the data shown in 
the right panels of Figs. 5 and 6 indicates that the most 
reliable description of the data is obtained by the ansatz 
albPftc 1 , which seems reasonable as the term c 1 In L/L 
decreases slowly and may therefore still be important at 
the available system sizes. 

A word of caution is in order here regarding the mean¬ 
ing of y 2 . This quantity is usually normalized by (gaus- 
sian) statistical errorbars attached to the data and should 
therefore follow the y 2 distribution. In particular, this 
implies that y 2 /ndf for a perfect fit approaches unity and 
can not be smaller unless the model “overfits” statistical 
noise. Here, the situation is strikingly different as our 
data do not bear statistical errorbars and y 2 does not 














5 




Figure 5. (Color online) Comparison of different fits over the 
range [L m i n ,10°] to the spin wave result for Si at J 2 = — 1 
in the Ji — J 2 Heisenberg model. The left panel displays the 
prefactor h of the logarithmic scaling term h ln(L), which all 
fits find to be very close to unity. The right panel displays 
the corresponding y 2 normalized by the number of degrees of 
freedom (ndf). Clearly the best fits with the lowest y 2 find 
h to be closest to 1. The artifacts around L m i n « 10 4 stem 
from a change of the grid on which we calculated Si which 
effectively introduces a higher weight for points in the denser 
region of the grid at smaller system size. The fit functions are 
coded according to the terms in equation (3.8). 


have any statistical meaning. In fact, for a perfect fit, y 2 
would then vanish, a situation we are very close to. The 
slow growth of the different fitting terms as well as the 
fact that the exact form of the subleading terms in the 
scaling below the logarithmic term remain unknown still 
gives rise to a small uncertainty of our fit results. 

Despite this, the different results for the investigated 
ansdtze consistently yield a logarithmic prefactor which 
is very close to (or evolves with growing system sizes into) 
l q = 1 for the Ji — J 2 Heisenberg model and l q = 1/2 for 
the XY model. This can be clearly seen in Figs. 5 and 6 
for ii, and for l q in Figs. 3 and 4 for different values of q 
(as well as different J 2 for the J\ — J 2 Heisenberg model). 

Our high-precision spin-wave results for a line subsys¬ 
tem are therefore in full agreement with the prediction 
Eq. (1.2) of a prefactor l q = no /2 reflecting the num¬ 
ber of Goldstone modes associated with the breaking of 
a continuous symmetry. 

From the structure of the eigenvalues of the correla¬ 
tion matrix Eq. (3.5) one can go further to interpret the 
additive logarithmic term in terms of the number of Gold- 
stone modes ng. Indeed, one can rewrite them as 




(3.9) 


Figure 6. (Color online) Comparison of different fits over 
the range [L m in, 7 • 10 4 ] to the spin wave result for Si in the 
XY model. The left panel displays the prefactor h of the 
logarithmic scaling term h ln(L), which all fits find to be very 
close to one half. The right panel displays the corresponding 
y 2 normalized by the number of degrees of freedom (ndf). 
Clearly the best fits with the lowest y 2 find h to be closest 
to 0.5. The fit functions are coded according to the terms in 
equation (3.8). 


with 


0(p, k y ) 


A[p , k y ) - B(p,k y ) 

A(p, ky) + B(p , ky) ’ 


(3.10) 


A and B being given by Eqs. (3.6) and (3.7), and the 
L modes p = —n + — 1) with j = 1,... ,L. It is 

straightforward to see that all 0 are non-singular 0(1) 
numbers, except at the singular points where Goldstone 
modes vanish. More precisely for SU(2) there are two 
contributions 


0 SU(2) (O,O) 


1 

0SU(2) ( 7 ^ tt) 



21Vtoaf, 


(3.H) 

where ttiaf is the thermodynamic limit (SU(2) broken) 
staggered magnetization, and one contribution for U(l) 


© U(1) (0,0) =* ^1? ^ 47Vm xy , (3.12) 

where m xy is the transverse order in the thermodynamic 
limit. Therefore all eigenvalues u p are 0(1) away from 
the Goldstone points where instead 

^Goldstone oc \/Z + constant. (3.13) 

Plugging this into the expression of the Renyi EEs 
Eq. (3.3), the L — ng modes with 0(1) eigenvalues will 
add up and contribute ~ L (the area law part) to S q and 
the reg terms will each contribute \ In L, \/q. 
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IV. DISCUSSIONS AND CONCLUSIONS 


We have investigated predictions from field theory that 
spontaneous breaking of a continuous symmetry leads to 
a logarithmic subleading scaling of the EEs S q indepen¬ 
dent on microscopic parameters and the Renyi index q , in 
the specific case of a periodic line subsystem embedded 
in a two-dimensional torus. Our results, obtained using 
two different methods (numerically exact QMC and spin 
wave theory), are in perfect agreement with the predic¬ 
tion that the prefactor of the logarithmic term is given 
by l q = no /2 by studying two models breaking SU(2) 
and U(l) symmetry respectively. 

Interestingly, we find that it is not necessary to study 
a bipartition of the system in two equal parts as cutting 
out a one dimensional subsystem is sufficient to capture 
the universal logarithmic correction. This is beneficial 
for both methods used in this work and we believe that 
other numerical and analytical techniques can profit from 
this finding in order to push calculations to larger system 
sizes, which are of tremendous importance for fitting the 
logarithmic term. Moreover, the spin-wave theory of the 
entanglement entropy of a line subsystem allows simpli¬ 
fied calculations where the contribution of each Gold- 
stone mode can be fully understood analytically in the 
modified (symmetry restored) spin-wave theory formal¬ 


ism. 

Reaching very large system sizes allowed us to capture 
higher order finite size corrections which demonstrates 
that it is very difficult to get a precise and size-converged 
estimate for the prefactor of the logarithmic correction l q 
using QMC simulations, restricted to linear sizes of a few 
tens of sites. 

Beyond this case of continuous symmetry-breaking 
phases, it would be interesting to investigate whether the 
line subsystem can also capture subdominant universal 
corrections associated with other types of phases, such 
as discrete symmetry-breaking or topological phases. In¬ 
deed in the latter case, a one-dimensional geometry, as 
used in Ref. 23, appears computationally more tractable 
(especially within QMC) than the usual topological en¬ 
tanglement entropy constructions 5,6,24 . 
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